Generalized Inverse Gaussian: Generalized Inverse Gaussian Distribution

Description

Density function, cumulative distribution function, quantile function and random number generation for the generalized inverse Gaussian distribution with parameter vector param. Utility routines are included for the derivative of the density function and to find suitable break points for use in determining the distribution function.

Usage

dgig(x, chi = 1, psi = 1, lambda = 1,
     param = c(chi, psi, lambda), KOmega = NULL)
pgig(q, chi = 1, psi = 1, lambda = 1,
     param = c(chi,psi,lambda), lower.tail = TRUE,
     ibfTol = .Machine$double.eps^(0.85), nmax = 200)
qgig(p, chi = 1, psi = 1, lambda = 1,
     param = c(chi, psi, lambda), lower.tail = TRUE,
     method = c("spline", "integrate"),
     nInterpol = 501, uniTol = 10^(-7),
     ibfTol = .Machine$double.eps^(0.85), nmax =200, ...)
rgig(n, chi = 1, psi = 1, lambda = 1,
     param = c(chi, psi, lambda))
rgig1(n, chi = 1, psi = 1, param = c(chi, psi))
ddgig(x, chi = 1, psi = 1, lambda = 1,
      param = c(chi, psi, lambda), KOmega = NULL)

Value

dgig gives the density, pgig gives the distribution function, qgig gives the quantile function, and rgig

generates random variates. rgig1 generates random variates in the special case where $\lambda = 1$.

ddgig gives the derivative of dgig.

Arguments

x, q: Vector of quantiles.
p: Vector of probabilities.
n: Number of observations to be generated.
chi: A shape parameter that by default holds a value of 1.
psi: Another shape parameter that is set to 1 by default.
lambda: Shape parameter of the GIG distribution. Common to all forms of parameterization. By default this is set to 1.
param: Parameter vector taking the form c(chi, psi, lambda) for rgig, or c(chi, psi) for rgig1.
method: Character. If "spline" quantiles are found from a spline approximation to the distribution function. If "integrate", the distribution function used is always obtained by integration.
lower.tail: Logical. If TRUE, probabilities are $P(X\leq x)$, otherwise as $P(X>x)$.
KOmega: Sets the value of the Bessel function in the density or derivative of the density. See Details.
ibfTol: Value of tolerance to be passed to incompleteBesselK by pgig.
nmax: Value of maximum order of the approximating series to be passed to incompleteBesselK by pgig.
nInterpol: The number of points used in qgig for cubic spline interpolation (see splinefun) of the distribution function.
uniTol: Value of tol in calls to uniroot. See uniroot.
...: Passes arguments to uniroot. See Details.

Author

David Scott d.scott@auckland.ac.nz, Richard Trendall, and Melanie Luen.

Details

The generalized inverse Gaussian distribution has density

$$f(x)=\frac{(\psi/\chi)^{\frac{\lambda}{2}}}% {2K_\lambda(\sqrt{\psi\chi})}x^{\lambda-1}% e^{-\frac{1}{2}\left(\chi x^{-1}+\psi x\right)}$$

for $x>0$, where $K_\lambda()$ is the modified Bessel function of the third kind with order $\lambda$.

The generalized inverse Gaussian distribution is investigated in detail in Jörgensen (1982).

Use gigChangePars to convert from the $(\delta,\gamma)$, $(\alpha,\beta)$, or $(\omega,\eta)$ parameterizations to the $(\chi,\psi)$, parameterization used above.

pgig calls the function incompleteBesselK from the package DistributionUtils to integrate the density function dgig. This can be expected to be accurate to about 13 decimal places on a 32-bit computer, often more accurate. The algorithm used is due to Slavinsky and Safouhi (2010).

Calculation of quantiles using qgig permits the use of two different methods. Both methods use uniroot to find the value of $x$ for which a given $q$ is equal $F(x)$ where $F$ denotes the cumulative distribution function. The difference is in how the numerical approximation to $F$ is obtained. The obvious and more accurate method is to calculate the value of $F(x)$ whenever it is required using a call to pghyp. This is what is done if the method is specified as "integrate". It is clear that the time required for this approach is roughly linear in the number of quantiles being calculated. A Q-Q plot of a large data set will clearly take some time. The alternative (and default) method is that for the major part of the distribution a spline approximation to $F(x)$ is calculated and quantiles found using uniroot with this approximation. For extreme values (for which the tail probability is less than $10^{-7}$), the integration method is still used even when the method specifed is "spline".

If accurate probabilities or quantiles are required, tolerances (intTol and uniTol) should be set to small values, say $10^{-10}$ or $10^{-12}$ with method = "integrate". Generally then accuracy might be expected to be at least $10^{-9}$. If the default values of the functions are used, accuracy can only be expected to be around $10^{-4}$. Note that on 32-bit systems .Machine$double.eps^0.25 = 0.0001220703 is a typical value.

Generalized inverse Gaussian observations are obtained via the algorithm of Dagpunar (1989).

References

Dagpunar, J.S. (1989). An easily implemented generalised inverse Gaussian generator. Commun. Statist. -Simula., 18, 703--710.

Jörgensen, B. (1982). Statistical Properties of the Generalized Inverse Gaussian Distribution. Lecture Notes in Statistics, Vol. 9, Springer-Verlag, New York.

Slevinsky, Richard M., and Safouhi, Hassan (2010) A recursive algorithm for the G transformation and accurate computation of incomplete Bessel functions. Appl. Numer. Math., In press.

Examples

Run this code

param <- c(2, 3, 1)
gigRange <- gigCalcRange(param = param, tol = 10^(-3))
par(mfrow = c(1, 2))
curve(dgig(x, param = param), from = gigRange[1], to = gigRange[2],
      n = 1000)
title("Density of the\n Generalized Inverse Gaussian")
curve(pgig(x, param = param), from = gigRange[1], to = gigRange[2],
      n = 1000)
title("Distribution Function of the\n Generalized Inverse Gaussian")
dataVector <- rgig(500, param = param)
curve(dgig(x, param = param), range(dataVector)[1], range(dataVector)[2],
      n = 500)
hist(dataVector, freq = FALSE, add = TRUE)
title("Density and Histogram\n of the Generalized Inverse Gaussian")
DistributionUtils::logHist(dataVector, main =
        "Log-Density and Log-Histogram\n of the Generalized Inverse Gaussian")
curve(log(dgig(x, param = param)), add = TRUE,
      range(dataVector)[1], range(dataVector)[2], n = 500)
par(mfrow = c(2, 1))
curve(dgig(x, param = param), from = gigRange[1], to = gigRange[2],
      n = 1000)
title("Density of the\n Generalized Inverse Gaussian")
curve(ddgig(x, param = param), from = gigRange[1], to = gigRange[2],
      n = 1000)
title("Derivative of the Density\n of the Generalized Inverse Gaussian")

Run the code above in your browser using DataLab